MoranI_0607
Ruby Xiong
7/6/2022
1 Moran’s I vs window size
stromaLabels = c("CD20pos Bcells", "CD4posTcell", "CD8posTcell", "CollagenCD45", "Stromal something", "Treg", "Macrophage")
epiLabels = levels(sc_HTAfull$label)[ !levels(sc_HTAfull$label) %in% stromaLabels]
sc_HTAfull$epiCell = ifelse(sc_HTAfull$label %in% stromaLabels, 0, 1)
sc_HTAfull[, levels(sc_HTAfull$label)] = model.matrix( ~ -1 + label, data=sc_HTAfull)
dat <- filter(sc_HTAfull, epiCell==1)
dat <- filter(dat, region=="region_001")
dat <- filter(dat, region=="region_001")
slides <- unique(dat$SlideID)[1:10]
dat <- filter(dat, SlideID %in% slides)
ADSSL <- c(2,3,5,7,9,1,4,6,8,10)
rs <- c(50, 100, 250, 500, 600)
markers <- c("TumorEpiMuc2","TumorStem","IEL","Goblet")
## plot function: Moran's I
moransI.df.transfer <- function(result.df){
result.df$marker1 <- rep(c("TumorEpiMuc2","TumorStem","IEL","Goblet"), each=4)
result.df$marker2 <- rep(c("TumorEpiMuc2","TumorStem","IEL","Goblet"), 4)
plot.df <- melt(result.df,variable.name = "radii", id.vars = c("marker1","marker2","slide","TumorType"),
value.name = "MoransI")
plot.df$radii <- as.numeric(levels(plot.df$radii))[plot.df$radii]
plot.df
}
plot.moransI <- function(result.df){
plot.df <- moransI.df.transfer(result.df)
p <- ggplot(plot.df, aes(x=radii, y=MoransI, group=slide)) +
geom_line(aes(color=TumorType), alpha=0.5) +
geom_point(shape=21, color="black", size=0.5, alpha=0.3)+ theme_bw()+
facet_grid(rows=vars(marker1), cols=vars(marker2), scales = "free_y")
p
}
# plot function for noramlized moran's I
mI.standard.transfer <- function(result.df){
result.df$marker <- c("TumorEpiMuc2","TumorStem","IEL","Goblet")
plot.df <- melt(result.df,variable.name = "radii", id.vars = c("marker","slide","TumorType"), value.name = "MoransI")
plot.df$radii <- as.numeric(levels(plot.df$radii))[plot.df$radii]
plot.df
}
plot.moransI.standard <- function(result.df){
plot.df <- mI.standard.transfer(result.df)
p <- ggplot(plot.df, aes(x=radii, y=MoransI, group=slide)) +
geom_line(aes(color=TumorType), alpha=0.5) +
geom_point(shape=21, color="black", size=0.5, alpha=0.3)+ theme_bw()+
facet_wrap(~marker)
p
}# plots_1 <- plots_2 <- list()
result1.all <- result2.all <-result3.all <-result4.all <-list()
for(j in ADSSL){
result1.df <- data.frame(matrix(NA, ncol=length(rs), nrow=16)) # moran's I
result2.df <- data.frame(matrix(NA, ncol=length(rs), nrow=4)) # Moran's I standardized
colnames(result1.df) <- colnames(result2.df) <- as.character(rs)
result1.df$slide <- result2.df$slide <- slides[j]
result3.list <- result4.list <- list()
for(i in 1:length(rs)){
result3.df <- data.frame(matrix(NA, ncol=4, nrow=sum(dat$SlideID==slides[j]))) # LISA
result4.df <- data.frame(matrix(NA, ncol=4, nrow=sum(dat$SlideID==slides[j]))) # LISA standardized
colnames(result3.df) <- colnames(result4.df) <-markers
result3.df$slide <- result4.df$slide <- slides[j]
dat_i <- filter(dat, SlideID == slides[j])
wlist = multiplexMoran::coords2listw(dat_i[,c("x", "y")], k = 500, maxDist = rs[i] )
moranMat = multiplexMoran::moran(select(dat_i, TumorEpiMuc2,TumorStem,IEL,Goblet),
listw = wlist, demean=TRUE)
result1.df[,i] <- as.vector(moranMat)
for (k in 1:length(markers)){
z_hat <- (moran.test(dat_i[, markers[k]], listw=wlist)$estimate)[1]
result2.df[k,i] <- z_hat
lisaVec = spdep::localmoran(dat_i[, markers[k]], listw = wlist)
result3.df[,k] <- lisaVec[,1]
result4.df[,k] <- lisaVec[,4]
}
result3.df$r <- result4.df$r <- rs[i]
result3.list[[i]] <- result3.df
result4.list[[i]] <- result4.df
}
result1.all[[j]] <- result1.df
result2.all[[j]] <- result2.df
result3.all[[j]] <- result3.list
result4.all[[j]] <- result4.list
# plots_1[[j]] <- plot.moransI(result1.df, name.slide = slides[j])
# plots_2[[j]] <- plot.moransI.standard(result2.df, name.slide = slides[j])
}
moran_result_list <- list(result1.all, result2.all, result3.all, result4.all)
save(moran_result_list, file="moran_0718.RData")setwd("/home/xiongj3/hist_fit_moran/atlasAnalysis")
load("moran_0718.RData")# for( i in 1:5 ){
# print(plots_1[[ADSSL[i]]])
# }
result1.all.df <- do.call("rbind", moran_result_list[[1]])
result1.all.df$TumorType <- rep(c("Adenoma", "SSL"), each=(nrow(result1.all.df)/2))
plot.moransI(result1.all.df)
# for( i in 6:10){
# print(plots_1[[ADSSL[i]]])
# }2 Standardized Moran’s I vs window size
# for( i in 1:5 ){
# print(plots_2[[ADSSL[i]]])
# }
result2.all.df <- do.call("rbind", moran_result_list[[2]])
result2.all.df$TumorType <- rep(c("Adenoma", "SSL"), each=(nrow(result2.all.df)/2))
plot.moransI.standard(result2.all.df)
# for( i in 6:10 ){
# print(plots_2[[ADSSL[i]]])
# }3 Number of cells v.s. Index
Effect of number of cells on value of Moran’s I
library(dplyr)
cell.num <- as.data.frame(table(dat$SlideID))
colnames(cell.num) <- c("slide", "cell.count")
plot.df <- moransI.df.transfer(result1.all.df)
cvi.plot.df <- left_join(plot.df, cell.num, by="slide")
ggplot(cvi.plot.df, aes(x=cell.count, y=MoransI))+
geom_point(aes(color=TumorType), alpha=0.5)+facet_wrap(~radii, nrow=1)+theme_bw()
Effect of number of cells on value of standardized Moran’s I
plot.df <- mI.standard.transfer(result2.all.df)
cvi.plot.df <- left_join(plot.df, cell.num, by="slide")
ggplot(cvi.plot.df, aes(x=cell.count, y=MoransI))+
geom_point(aes(color=TumorType), alpha=0.5)+facet_wrap(~radii, nrow=1)+theme_bw()
4 LISA
LISA is taken \(\sqrt[\leftroot{-1}\uproot{2}\scriptstyle 5]{LISA}\qquad\) for illustration purposes
result3.list <- moran_result_list[[3]]
adssl <- rep(c("Adenoma", "SSL"), each=5)
library(cowplot)
cube.root <- function(number){
sign <- (-1)^(1*(number<0))
sign*abs(number)^(1/5)
}
plot.lisa <- function(plot.lisaMat, root.trans=TRUE){
plot.lisa.df <- melt(plot.lisaMat,variable.name = "Markers", id.vars = c("slide","r","x","y"),
value.name = "LISA")
if(root.trans==TRUE){
plot.lisa.df$LISA <- cube.root(plot.lisa.df$LISA)
}
ggplot(data=plot.lisa.df, aes(x=x, y=y, color=LISA)) +
geom_point(size=0.05, alpha=0.3)+
scale_color_gradient2()+
theme_dark()+theme(legend.position="top",legend.text = element_text(angle = 50))+
facet_grid(rows=vars(r), cols=vars(Markers), scales = "free_y")
}
lisa.df.print <- function(i){
result3.df <- result3.list[[i]]
plot.lisaMat <- do.call("rbind", result3.df)
plot.lisaMat$x <- dat[dat$SlideID==plot.lisaMat$slide[1],"x"]
plot.lisaMat$y <- dat[dat$SlideID==plot.lisaMat$slide[1],"y"]
p1 <- plot.lisa(plot.lisaMat)+ggtitle(paste(slides[i], adssl[i]))
dat.i <- dat[dat$SlideID==plot.lisaMat$slide[1],]
dat.i <- subset(dat.i, dat.i$label %in% markers)
dat.i$label <- factor(dat.i$label, levels = markers)
p2 <- ggplot(dat.i, aes(x=x, y=y))+geom_point(size=0.05)+
theme_bw()+facet_wrap(~label, scales = "fixed", nrow=1)
pp <- plot_grid(p1, p2, nrow = 2, rel_heights = c(4/5, 1/5))
print(pp)
}4.1 HTA11_10167_0000_01_01
lisa.df.print(1)
4.2 HTA11_10623_0000_01_01
lisa.df.print(2)
4.3 HTA11_10711_0000_01_01
lisa.df.print(3)
4.4 HTA11_4255_0000_02_02
lisa.df.print(4)
4.5 HTA11_6298_0000_04A_03
lisa.df.print(5)